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^ ■ Abstract 
o . 

. The Monte Carlo wave function method or the quantum trajectory /jump approach is a powerful 

'sj" i tool to study dissipative dynamics governed by the Markovian master equation, in particular for 

high-dimensional systems and when it is difficult to simulate directly. In this paper, we extend 
Q_i' this method to the non-Markovian case described by the generalized Lindblad master equation. 

^ , Two examples to illustrate the method are presented and discussed. The results show that the 

a 

^ I method can correctly reproduce the dissipative dynamics for the system. The difference between 

this method and the traditional Markovian jump approach and the computational efficiency of this 
^ . method are also discussed. 
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I. INTRODUCTION 



Since the pioneering work of Albert Einstein, who explained the phenomenon of dissipa- 
tion and Brownian motion in his annus mirabilis of 1905 by use of statistical methods, a rich 
variety of methods to tackle quantum fluctuations and quantum dissipation in open systems 
has been proposed l|, |2|. Among them the quantum master equation (QME) approach and 
the quantum Langevin description (QLE)[3] are two of powerful functional integral tech- 
niques for the study of time evolution of open quantum systems. The quantum master 
equation can be divided into two categories: Markovian and non-Markovian. The Marko- 
vian master equation {4] (especially in the Lindbald form) can be derived with the weak 
coupling limit (or the Born approximation) and the Markovian approximation. It can be 
solved analytically!^ for some special cases, but for most cases we have to solve and simu- 



late it numerically by the Monte Carlo wave function method or quantum trajectory /jump 

. This method is very effective for qubit systems even with large 



approach 0, Q, [s, 
number of qubits, say n = 24[ 

However, the dynamics of an open system is not always Markovian. Strong system- 
environment couplings, correlation and entanglement in the initial state and structured 
reservoirs may lead the dynamics far from Markovian. Many methods have been proposed 
to describe the non-Markovian process, including the Lindblad equation with time dependent 
decay rates [3], generalized Lindblad equation [13}] obtained from the correlated projection 

, [isl, phenomenological memory kernel master equation [igI. \v\ 



and the post-Markovian master equation 



18, 



19 



20]. The first two methods are local in time 



while the last two involve an integral of time. For the first method, the only difference from 
the Markovian master equation is that the decay rates in the equation are time-dependent. 
These decay rates may take not only positive values but also negative ones. When decay 
rates are positive, the Markovian Monte Carlo wave function method can directly be used. 
However, the mehtod is not available when the decay rates are negative. This problem was 



solved in Ref. 



2l\ by introducing reversed jumps. 



The generalized Lindblad master equation can well describe the dynamics of an open sys- 
tem beyond the Markovian limit, especially it is very effective for an environment composed 



of spins 



22,123, 



24l | and structured reservoirs 251]. However the extension of the Monte Carlo 



simulation to this equation remains untouched. In this paper, we will explore the unravel- 
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ing and quantum trajectory approach for the generahzed Lindblad equation. The structure 
of this paper is organized as follows. In SecjlT] we briefly review the generalized Lindblad 
equation. In Sec JIIII we give the unraveling of this equation and generalize the Monte Carlo 
method to this equation. Two examples are presented in Sec JIVI Finally, we conclude our 
results in SecJVl 



II. GENERALIZED LINDBLAD MASTER EQUATION 

The equation that governs the dynamics of an ope n q uantum system can be derived 
by means of the projection superoperator technique 13, llJ]- The form (Markovian or non- 
Markovian) of the master equation crucially depends on the approximation used in the 
derivation, reflecting in the projection superoperator chosen. When we project the total 
system state into a tensor product, we can obtain the Markovian master equation, whereas 
a non-Markovian master equation can be obtained when we use a correlated projection. 
The following is the master equation derived by this method and it is called the generalized 



Lindblad master equation [l3| 



where Hm are Hermitian operators and R^n arbitrary system operators depending on 
the form of system-environment interactions. If we have only a single component ps = Pi, 
this equation obviously reduces to the ordinary Markovian master equation. In this paper 
we will focus on the case where we have at least two components. The state of the reduced 
system in this case is ps = J2mPm, we remind that Trp^ < 1- 



III. QUANTUM JUMP 

For clarity, we define the jump operators W^^ = and non-jump operators 

— ^ ~ iTirndt, whcrc the non-Hermitian effective Hamiltonian is given by Tim = 
Hm — ^iJ^nxR-nlnRnm- There are two subscripts and one superscript for the operator W^^. 
The first subscript m denotes the index of component where the system in, while the second 
subscript n denotes the index of component for the operation acting on; the superscript A 
represents the jump mode. Initially we assume that each operator Pm(^o) can be written 
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as Pm(to) = \'ipmi'to)){4'm(to)\, wheie \ipm(to)) is a non-normalized wave function. After an 
infinitesimal time dt, it evolves into the following state 



Pmito + dt) = J2 \^L){^L\dpL + l^lJi^lJdplm, (2) 
where the new states are defined by 



and 



with probabilities 



|,0 \ _ V^^mmlV^m(t)) 



dpL = ^{Mto)\W^nWL\Mto))dt, 

Pm 

dplm = —{^mitoW^mW^J^/jM), (5) 
Pm 

respectively. In Eqs.dl]) and (jl]), 

nX 

is the weight for the component pm that satisfies 

Pm = TTpm(t + dt). (7) 

Note that the jumps for pm depend on the other components pn, {n ^ m) of the reduced 
density matrix p. This makes our method different from the traditional quantum jump 
method. 

We can prove this unraveling by taking the jump and non-jump states ([3]), (jl]) and the 
probabilities dHD into Eq.®, 

Pm{t + dt) 

= E w:;,Mt)){i^n{t)\w±dt + w'^J^m{t)){i^m{t)\w^J^. (8) 

Simple algebra shows that in the limit dt 0, Eq([H]) reveals Eq.([T]). The evolution governed 
by Eq.([T]) can be simulated numerically by the so-called Monte Carlo wave function approach 
according to the unraveling given above. We start the time evolution from the state p(to) = 



J2mPm(to) = Em I V'm(^o)) (V'm(io) I , wheie Pm("^ = 1,2,3, ...) are the components for p. At 
time to + dt, where dt is much smaller than the timescale relevant for the evolution of the 
density matrix, a random number e which is randomly distributed in the unit interval [0, 1] 
is used to determine the jump. Note that all the components are controlled by this random 
number. For each component \ipm), if < e < dpj^i, it jumps to iV'mi), if c^Pmi < e < 
dPmi + dPmiy it jumps to iV'^i), and so on. These jumps are all operated on the component 
Pi; if J2xdp^i < e < T.xdPmi + ^Pm2' it jumps to the component 2, namely l^/'^a)- Jumps 
to the other components can be established in a similar way. If e > EnA dPmm ^ non-jump 
takes place and the state ends up in iV'mm)- This operation is acted on the component pm 
itself. We define a generalized jump superoperator Wj, which denotes all jumps for all the 
components controlled by this random number. We repeat this process as many times as 
n = At/dt for all the components, where At is the total evolution time. We call this single 
evolution a generalized quantum trajectory. This trajectory contains all the components of 
the density matrix. Given an operator A, we can write its mean value {A){t) = Tr{Ap{t)) 
as an average over J\f trajectories as 

{A){t) = jim^^5:(^„,,.(t)|A|^„,,(t)). (9) 

j=l m 

IV. APPLICATION 



In this section, we use the model and the generalized master equation given in Refs. 25| 

and [23] as two examples to illustrate our method. First consider a two-state system coupled 

to an environment. The environment consists of a large number of energy levels which are 

arranged into two energy bands with the same energy spacing(see FigH]). The lower energy 

band contains A^^i levels while the upper one N2 levels. This model can be understood 

as a "many level" environment or "container", of which only the relevant parts of the 

spectrum enter the model. For details of this model, we refer the reader to H]- The 

total Hamiltonian for a qubit coupled to such an environment in Schrodinger picture is 

H = Ho + V[25] with(we set h = I) 

1 S€ 5^ 

Ho = + + — n2)|n2)(n2|, 

V = c(ni, n2)o-+|ni)(n2| +H.c., 
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where the index ni denotes the levels of lower energy band and n2 denotes the levels of upper 
band, and are Pauli operators. A is the overall strength of the interaction, c{ni,n2) 
are coupling constants, they are independent of each other and are identically distributed, 
satisfying 

(c(?2i,n2)) = 0, 
{c{ni,n2)c{n[,n'^)) = 0, 
{c{ni,n2)c*{n[,n2)) = 5„i,„;5„2,n^- 

According to Hq, one can transform the problem into the interaction picture and, with the 
help of projection superoperator technique, obtain the non-Markovian evolution equation as 

j/s\t) = li<y^pf{t)a- - |{a+a-,pW(t)}, 

j/ht) = l2cr-pf{t)a^ - |{a-a+,p(^)(t)}, (10) 

where 

li = — j-^ — (^ = 1,2). 

With definitions of IIi = and 112 = Z]n2 1^2) (^^2 1, 111 + 112 = Ie^ the two 

non-normalized density matrixes can be obtained by p^^ = TTE^^iPT)-,^ = 1,2, where pt 
is the total density matrix for the system and environment. The reduced density matrix 
for the system is then given by p = pg^ + p^g\ We note that in Eq. (|TOl) . there are no 
environment operators other than the two (c-number) parameters 71,72- The initial state 
of the environment is taken into account by means of the distribution of initial Ps^p^^ its 
effect on the system dynamics was plotted in Figsl2] and [31 This equation can be written 
in the form of Eq.([T]) by setting Hi = 0, i?ii = R22 = 0, i?i2 = ^Ticr"'', and R21 = ^/^y2<^~ ■ 
In this model, there is only one jump operator for each component, i.e. = V^*^^ 

and W21 = \f^(y~ , and non-jump operators = I — iHmdt with Hi = — |72(T^cr^ and 

^2 = -^Jicr'cr+. 

We consider two types of initial condition in the following simulation. First, only the 

(2) 

lower band of the environment is populated, i.e. pg =0. Under this condition, the reduced 
system can be solved analytically. Another case is, the two bands of the environment are 
all populated. With this initial condition, we solve the master equation numerically. In 
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both cases, we choose an initial states |0(O)) = |e) and |0(O)) = |(|e) + \g)) for the system, 
where |e) and \g) denote the excited state and ground state, respectively. We compare the 
analytic solution and the numerical simulation (solve the equation by Runge-Kutta method) 
to the results obtained from the quantum jump/trajectory approach in Figsj2] and [3l The 
trajectory number in this quantum jump approach is A/" = 400. We can see from the figures 
that the quantum trajectory approach correctly reproduces the system evolution. The errors 
are sufficiently small, although we choose a small number of trajectories, showing that this 
method is efficient. 



Another example is a qubit coupled to a spin bath[23|. The full system consists of a 



central spin interacting with a bath of N spins. Such a system can be described by 

u ^ 

H = —a^ + akB ■ ffk, (11) 
k=i 

where a denotes the Pauli matrix for the central spin that is the system we are interested in, 
and (Tk stands for the k-th spin in the bath. After defining an unperturbed part Hq = + 
2azKz, where = ^Y,k=i'^kO'z, the Hamiltonian can be transformed into the interaction 
picture. Assuming the parameters are real and time independent, the master equation reads 

where pm = T^i^B{PT^m), Pt is the density matrix for the total system(the central spin 
plus the bath), 11^ is a projection superoperator that projects the z-component of the 
bath angular momentum into an eigenvector with eigenvalue m. We take = 2 as an 
example, then the density matrix of the central spin has three components, denoted by 
pi,P0)P-i; respectively. Each component has two jump operators which act on the other 
two components, and a non-jump operator, which acts on itself. The comparison between 
directly numerical simulations (by Runge-Kutta) and quantum trajectory method is shown 
in FigJH Here the trajectory number is chosen to be A/" = 4000. We can find that as the 
number of jump operators and components increases, the number of quantum trajectory, 
with that we can obtain a correct result, increases. 



V. CONCLUSION AND DISCUSSION 



In this paper, we have developed an efficient unraveling for the generalized Lindblad 
master equation. Based on this unraveling, a generalized Monte Carlo wave function method 
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is presented. It is worth addressing that in this Monte Carlo wave function method, we need 
only to store M non-normalized wave function, i.e. M length-A^ vectors (M denotes the 
number of the components for the reduced density matrix and stands for the dimension of 
the Hilbert space) instead of the density operator, which are M N x N matrices, hence this 
method saves the computer time and space. The difference between the ordinary quantum 
jump method and the present one is that the latter describes a non-Markovian dynamics. 
In addition, the point that each component pi of the density matrix is non-normalized, and 
jumps along the component depend on the component other than pi is also different. 
By successfully simulating the coupling among those components, this method can simulate 
the non-Markovian dynamics efficiently. Further examination shows that the computational 
complexity increases with the number of the components. The increased complexity due to 
the increase of the components and jump operators can by analyzed as follows. Assume the 
jump operators and the number of jump operators are restricted to be the same for each 
component, the possible jump mode for g = (pi,p2, ■ ■ ■), or the number of the generalized 
jump superoperators W is 

A = M(J-1) + 1. (13) 

Here J is the number of jump operators for each component (including the non-jump oper- 
ator). The role that A plays is similar to the number of the jump operators in the ordinary 
Markovian master equation. It is well known that one downside of the quantum jump ap- 
proach is the complexity growth as the jump operators proliferate. From Eq.( fT3|) . we can 
find that this downside still exists in the presented method. Still, our method is effective 
when one simulates the decoherence governed by the non-Markovian master equation, as 
well as for a system with Hilbert space of high dimension. 
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FIG. 1: A two-state system coupled to an environment consisting of two energy bands with a finite 
number of levels. 
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FIG. 2: (Color online) Comparisons of the analytic solution for Eq. (jlOp to the results given by 
quantum trajectory approach. The initial state of the system is chosen |0(O)) = |e) in the top figure 
while 10(0)) = "^(|e) + l^)) in the bottom one. Initially only the lower band of the environment is 
populated. The other parameters chosen are 71 = 72 = 1- The time t is plotted in units of 1/h. 
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FIG. 3: (Color online) Comparisons of the numerical solution for Eq. ijlOp to the results given by 
quantum trajectory approach. The initial state of the system is chosen |</>(0)) = |e) in the top 
figure while |</'(0)) = + l^)) in the bottom one. Initially the two bands of the environment 

are populated. Two parameters are 71 = 72 = 1- The time t is plotted in units of 1/h. Note that 
in the bottom figure we plot the off-diagonal element peg of the reduced system. 
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FIG. 4: (Color online) Comparisons of the numerical solution for Eq. ()12p to the results given by 
quantum trajectory approach. The initial state of the system is chosen pi = po = p-i = 
All parameters in the equation are set to be equal. The time t is plotted in units of 1/h. 
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